##### Replication script for "Who Counts? Non-Citizen Residents, Spatial Sorting, and Malapportionment"

### Packages
library(openxlsx)
library(lmtest)
library(sandwich)
library(texreg)
library(tidyverse)

### Data

all = read.xlsx("PATH.../mal.xlsx")

pr = read.xlsx("PATH.../malpr.xlsx")

### Reference Group Interaction Term

all$denom = factor(all$denom, , levels = c("Population", "Citizens"))
levels(all$denom)
pr$denom = factor(pr$denom, , levels = c("Population", "Citizens"))

### Estimation

# Table 1

m00 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter, all, weights = basew)
m00se <- coeftest(m00, vcovCL(m00, ~ID))[,2]
m00p <- coeftest(m00, vcovCL(m00, ~ID))[,4]

m1 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter, all)
m1se <- coeftest(m1, vcovCL(m1, ~ID))[,2]
m1p <- coeftest(m1, vcovCL(m1, ~ID))[,4]

m2 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter, pr, weights = basew)
m2se <- coeftest(m2, vcovCL(m2, ~ID))[,2]
m2p <- coeftest(m2, vcovCL(m2, ~ID))[,4]

texreg(list(m00,m1,m2), override.se = list(m00se,m1se,m2se),
       override.pvalues = list(m00p,m1p,m2p), digits = 4,
       omit.coef = "(ID)|(ctry)|(Year)",#stars = c(0.001, 0.01, 0.05, 0.1),
       custom.header = list("With Weights" = 1, "Without Weights" = 2, "PR Only" = 3),
       custom.coef.map = list("lmal" = "Lag RRI", "denomCitizens" = "Reference Citizens",
                              "excl" = "Share Non-citizens", "excl:denomCitizens" = "Share Non-citizens * Ref. Citizens",
                              "minmag" = "Below Min. Rep. Threshold","lvoter" = "ln Number of Voters"),
       booktabs = T, dcolumn = T, use.packages = F,custom.gof.rows = list(
         "Country FE" = c(rep("Yes",3)),
         "Year FE" = c(rep("Yes",3))),
       caption = "The Effect of Non-Citizen Residents and the Apportionment Mechanism on the Relative Representation Index (RRI)",
       label = "est1", custom.note = "Cluster-Robust Standard Errors by District. %stars") 

# Figure 2

# Like the specification in the first model in table 1 but alternative formula  to get the point estimates directly

mplot = lm(mal~lmal+denom/excl+factor(ctry)+factor(Year)+minmag+lvoter, all, weights = basew)

plotdat <- data.frame(coeftest(mplot, vcovCL(mplot, ~ID))[c("denomPopulation:excl","denomCitizens:excl"),c(1:2)])

plotdat <- plotdat %>% rownames_to_column(., var = "group") %>%
  mutate(conf.low = Estimate-Std..Error*1.96,conf.high = Estimate+Std..Error*1.96,
         group = ifelse(group == "denomPopulation:excl", "Population","Citizens"))

ggplot(plotdat, aes(group, Estimate))  + theme_bw()+
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .05)+ geom_point() + geom_hline(yintercept = 0)+
  labs(x = "Apportionment Mechanism", y = "Marginal Effect of the Share of Non-Citizen Residents")+
  theme(axis.text=element_text(size=12),axis.title=element_text(size=14,face="bold")) 

# Table A5

m4 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter, all, weights = basew)
m4se <- coeftest(m4, vcovCL(m4, ~ID))[,2]
m4p <- coeftest(m4, vcovCL(m4, ~ID))[,4]

m5 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter, all)
m5se <- coeftest(m5, vcovCL(m5, ~ID))[,2]
m5p <- coeftest(m5, vcovCL(m5, ~ID))[,4]

m6 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter, pr, weights = basew)
m6se <- coeftest(m6, vcovCL(m6, ~ID))[,2]
m6p <- coeftest(m6, vcovCL(m6, ~ID))[,4]

texreg(list(m4,m5,m6), override.se = list(m4se,m5se,m6se),
       override.pvalues = list(m4p,m5p,m6p),  digits = 4,
       omit.coef = "(ID)|(ctry)|(Year)",#stars = c(0.001, 0.01, 0.05, 0.1),
       custom.header = list("With Weights" = 1, "Without Weights" = 2, "PR Only" = 3),
       custom.coef.map = list("lmal" = "Lag RRI", "denomCitizens" = "Reference Citizens",
                              "excl" = "Share Non-citizens", "excl:denomCitizens" = "Share Non-citizens * Ref. Citizens",
                              "minmag" = "Below Min. Rep. Threshold","lvoter" = "ln Number of Voters"),
       booktabs = T, dcolumn = T, use.packages = F,custom.gof.rows = list(
         "Country FE" = c(rep("Yes",3)),
         "Year FE" = c(rep("Yes",3)),
         "District FE" = c(rep("No",0),rep("Yes",3))),
       caption = "The Effect of Non-Citizen Residents and the Apportionment Mechanism on the Relative Representation Index (RRI): With District FE",
       label = "est2", custom.note = "Cluster-Robust Standard Errors by District. %stars") 

# Table A6


m00 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter+app, all, weights = basew)
m00se <- coeftest(m00, vcovCL(m00, ~ID))[,2]
m00p <- coeftest(m00, vcovCL(m00, ~ID))[,4]

m1 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter+app, all)
m1se <- coeftest(m1, vcovCL(m1, ~ID))[,2]
m1p <- coeftest(m1, vcovCL(m1, ~ID))[,4]

m2 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+lvoter+app, pr, weights = basew)
m2se <- coeftest(m2, vcovCL(m2, ~ID))[,2]
m2p <- coeftest(m2, vcovCL(m2, ~ID))[,4]

m4 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter+app, all, weights = basew)
m4se <- coeftest(m4, vcovCL(m4, ~ID))[,2]
m4p <- coeftest(m4, vcovCL(m4, ~ID))[,4]

m5 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter+app, all)
m5se <- coeftest(m5, vcovCL(m5, ~ID))[,2]
m5p <- coeftest(m5, vcovCL(m5, ~ID))[,4]

m6 <- lm(mal~lmal+excl*denom+factor(ctry)+factor(Year)+minmag+ID+lvoter+app, pr, weights = basew)
m6se <- coeftest(m6, vcovCL(m6, ~ID))[,2]
m6p <- coeftest(m6, vcovCL(m6, ~ID))[,4]

texreg(list(m00,m1,m2,m4,m5,m6), override.se = list(m00se,m1se,m2se,m4se,m5se,m6se),
       override.pvalues = list(m00p,m1p,m2p,m4p,m5p,m6p),  digits = 4,
       omit.coef = "(ID)|(ctry)|(Year)",#stars = c(0.001, 0.01, 0.05, 0.1),
       custom.header = list("With Weights" = 1, "Without Weights" = 2, "PR Only" = 3, "With Weights" = 4, "Without Weights" = 5, "PR Only" = 6),
       custom.coef.map = list("lmal" = "Lag RRI", "denomCitizens" = "Reference Citizens",
                              "excl" = "Share Non-citizens", "excl:denomCitizens" = "Share Non-citizens * Ref. Citizens",
                              "minmag" = "Below Min. Rep. Threshold","lvoter" = "ln Number of Voters", "app" = "Apportionment"),
       booktabs = T, dcolumn = T, use.packages = F,custom.gof.rows = list(
         "Country FE" = c(rep("Yes",6)),
         "Year FE" = c(rep("Yes",6)),
         "District FE" = c(rep("No",3),rep("Yes",3))),
       caption = "The Effect of Non-Citizen Residents and the Apportionment Mechanism on the Relative Representation Index (RRI): With Apportionment Control",
       label = "appor", custom.note = "Cluster-Robust Standard Errors by District. %stars") 
